clc;
clear;

domain_type = "tri";

switch domain_type
    case "1D"
        gauss_points = [-0.9681602395076260898355762, -0.8360311073266357942994298, -0.6133714327005903973087020, -0.3242534234038089290385380, 0, 0.3242534234038089290385380, 0.6133714327005903973087020, 0.8360311073266357942994298, 0.9681602395076260898355762];
        gauss_weights = [0.0812743883615744119718922, 0.1806481606948574040584720, 0.2606106964029354623187429, 0.3123470770400028400686304, 0.3302393550012597631645251, 0.3123470770400028400686304, 0.2606106964029354623187429, 0.1806481606948574040584720, 0.0812743883615744119718922];
    case "tri"
        gauss_points = [(1-sqrt(3/5))/2, (1+sqrt(3/5))*(1-sqrt(3/5))/4;
            (1-sqrt(3/5))/2, (1+sqrt(3/5))/4;
            (1-sqrt(3/5))/2, (1+sqrt(3/5))*(1+sqrt(3/5))/4;
            1/2, (1-sqrt(3/5))/4;
            1/2, 1/4;
            1/2, (1+sqrt(3/5))/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))*(1-sqrt(3/5))/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))*(1+sqrt(3/5))/4];
        gauss_weights = [100/324*(1+sqrt(3/5))/8, 5/81*(1+sqrt(3/5)), 100/324*(1+sqrt(3/5))/8, 5/81, 8/81, 5/81, 100/324*(1-sqrt(3/5))/8, 5/81*(1-sqrt(3/5)), 100/324*(1-sqrt(3/5))/8];
    case "rec"
        gauss_points = [-sqrt(3/5), -sqrt(3/5);
            -sqrt(3/5), 0;
            -sqrt(3/5), sqrt(3/5);
            0, -sqrt(3/5);
            0, 0;
            0, sqrt(3/5);
            sqrt(3/5), -sqrt(3/5);
            sqrt(3/5), 0;
            sqrt(3/5), sqrt(3/5)];
        gauss_weights =[25/81, 40/81, 25/81, 40/81, 64/81, 40/81, 25/81, 40/81, 25/81];
end

fprintf("domain_type: \n\t%s\n\n", domain_type);
fprintf("gauss_points:\n");
vpa(gauss_points)
fprintf("gauss_weights:\n");
vpa(gauss_weights)
